% ---------------------------------------------------------------------------
%   2-D, MULTIGROUP DIFFUSION-BASED, NONLINEAR RESPONSE MATRIX SOLVER - TEST!
% ---------------------------------------------------------------------------
%  j. roberts, 01/23/2010
%   This MAIN file currently only drives a power iteration like scheme.  As
%   soon as possible, the nonlinear solver will be used.  This file depends
%   on:
%       redblack2D.m      - iterative solver
%       resp2dexample.m   - response function generatory (i.e. a diff. code)
%       twoDcoefMGresp.m  - coefficients for diff solve
%       cresp.m           - current responses
% ---------------------------------------------------------------------------

% 2 --> 0.60307
% 25 --> 0.04%%%
% 50 -->

clear; clc;
format long
global order numg numel

%----------------------- GET PROBLEM DATA -------------------------
for or = 0:7
order = or; numg = 2; numel = 2;
j = ones(1,numel*4*(order+1)*numg); % a uniform guess of coefficients
j = j / sqrt(j*j');
if numel == 1
    kk= 0.03888561054922; % reference keff
else
    kk= 1.13977063757579;
end
x = [j kk]';
tol = [1.e-8 1.e-8]; % tolerance for norms
[x,normk,it,nrm] = redblack2D(x);
disp([num2str(or),'  ', num2str( x(end) ), '   ', ...
        num2str(100*(kk-x(end))/x(end))])

end
return
%---------------------- NEWTON SOLVER -----------------------------
% uses the standard newton method w/ a fd-approximated jacobian
disp('----NEWTON----')
gm = 0;
tic
x3          = [x' 1.0]'; 
xx3(:,1)    = x3;
z           = respfct2Dexample(x3); % initial residual
itmx        = 20; 
it          = 1;
zz(1)       = norm(z);
nrmcn(1)    = 1;
while ( zz(it) > tol(1)*zz(1)+tol(2) && it < itmx )
    fp          = jacob('respfct2Dexample',x3,z);
    s           = -fp\z;
    x3          = x3 + s;
    it          = it+1; 
    xx3(:,it)   = x3; % Keep all x's to estimate rho
    z           = respfct2Dexample(x3);
    zz(it)      = norm(z);
    x3(end-1)
    disp(zz(it))
end
t3 = toc;
plotcn = zz/zz(1); % relative nonlinear residual
%p3 = convest(xx3,xsol,0); % est. rate of conv. w/ RB soln. as ref
return
%--- OUTPUT
semilogy(0:itrb-1,plotrb,'k',0:itgm-1,plotgm,'b--',...
         0:it-1,plotcn,'r-.')
axis([0,max(itgm,max(itrb,it)), ...
    min(min(plotcn(end),plotgm(end)),plotrb(end)), 1.01]);     
xlabel('Iterations')
ylabel('Relative Nonlinear Residuals')
legend('red-black','newt-gmres','newton')

  disp(' ')
  disp(' *** final results ***')
  disp(' ')
  disp('         |    red-black    |    newt-gmres    |    classic      |   ')
fprintf('      it |      %3i        |      %3i         |     %3i         | \n', ...
    itrb-1,itgm-1,it-1)
fprintf('    keff |%16.13f |%16.13f  |%16.13f | \n', x1(end),x2(end-1),x3(end-1))

n2=norm( feval('respfctAUG',x2) );
n3=norm( feval('respfctAUG',x3) );
fprintf('||F(x)|| |    %4.3e   |    %4.3e    |    %4.3e   | \n',n1,n2,n3)
fprintf('       p |      %4.3f      |      %4.3f       |      %4.3f      | \n',p1,p2,p3)
fprintf('    time |      %4.3f      |      %4.3f       |      %4.3f      | \n',t1,t2,t3)
